% y代表原序列
[y,Fs]=audioread('test.wav');
y=y(:,1);
T=1/Fs;
t=(1:length(y))*T;
subplot(3,2,1)
plot(t,y)
ylim([-0.6 0.6])
xlabel("时间/s");
ylabel("幅度");
title("原始语音信号时域图");
Y=fft(y);
subplot(3,2,2)
plot((1:length(Y))*2/length(Y),abs(Y))
xlabel("频率");
ylabel("幅度");
title("原始语音信号频域图");
ylim([0 600])
% 加入高斯白噪声，即加噪后的序列y1
y1=y+0.2*rand(size(y))-0.1;
%sound(y1,Fs)
% audiowrite("加噪.wav",y1,Fs);
subplot(3,2,3)
plot(t,y1)
ylim([-0.6 0.6])
xlabel("时间/s");
ylabel("幅度");
title("加噪音语音信号时域图");
Y1=fft(y1);
subplot(3,2,4)
plot((1:length(Y1))*2/length(Y1),abs(Y1))
ylim([0 600])
xlabel("频率");
ylabel("幅度");
title("加噪音语音信号频域图");
% 降噪处理
% 维纳滤波去噪

% 计算信号和噪声的功率谱密度
Pyy = abs(Y).^2;
Pnn = abs(Y1 - Y).^2;

% 计算维纳滤波器的频率响应
H_wiener = Pyy./ (Pyy + Pnn);

% 对加噪信号进行维纳滤波
Y2_wiener = Y1.* H_wiener;
% sound(y2_wiener, Fs);
% 逆傅里叶变换得到时域去噪信号
y2_wiener = ifft(Y2_wiener);

% 去噪音语音信号时域图
subplot(3,2,5);
plot(t,y2_wiener);
ylim([-0.6 0.6]);
xlabel("时间/s");
ylabel("幅度");
title("维纳滤波去噪音语音信号时域图");


% 去噪音语音信号频域图
subplot(3,2,6);
plot((1:length(y2_wiener))*2/length(y2_wiener),abs(Y2_wiener));
ylim([0 600]);
xlabel("频率");
ylabel("幅度");
title("维纳滤波去噪音语音信号频域图");


% 计算滤波前的信噪比
% 原始信号功率

psignal=P(y);
% 噪声功率
pnoise=P(y1-y);
% 信噪比
SNR_Original=10*log10(psignal/pnoise);
disp(SNR_Original)

% 计算滤波后的信噪比
% 降噪后的功率
psignal_after=P(y);
% 噪声功率
pnoise=P(y2_wiener-y);
SNR_After=10*log10(psignal_after/pnoise);
disp(SNR_After)




